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JiB^RACT 

The  method  of  front  tracking  is  applied  to  problems  involving  curved 
detonation  fronts,  surface  instabilities  and  two-dimensional  Riemann  prob- 
lems.  The  detonation  problems  include  detonation  fronts  with  and  without 
cylindrical  symmetry;  comparisons  with  one-dimensional  models  are  made. 
The  analysis  of  interface  instabilities  focuses  on  the  compressible  Rayleigh- 
Taylor  instability  of  a  supersonic  accelerated  contact  discontinuity  between 
two  gases  and  the  propagation  of  a  supersonic  slab  jet.   Theoretical  notions 
for  an  S  matrix  theory  for  general  multi-dimensional  hyperbolic  conservation 
laws  and  the  numerical  implementation  of  computer  programs  which  solve 
certain  two-dimensional  Riemann  problems  are  also  discussed. 
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1.  Introdaction 


bangiaa.. 
•?ii£Jsb  /. 


Systems  of  non-linear  conservation  laws  in  n  space  dimensions 

a,  +  ?F(x,  u)  =  0  (1.1) 

where  o  =  u(x,  r)  and  F  is  a  smooth  function  of  the  state  a  e  R'  and  the  position  x  e  R"  into 
R^,  are  often  used  as  first  order  approximations  for  many  natural  phenomena.  Equations  of 
this  type  occur  in  models  in  which  external  forces  and  higher  order  effects  such  as  viscosity 
and  heat  conduction  are  neglected. 

We  are  primarily  concerned  with  the  case  where  system  (1.1)  consists  of  the  Euler 

equations  for  a  compressible,  inviscid,  non-heat  conducting  gas.  In  this  case: 

/  \ 

m 

mem 


D  = 


m     and  F(o)  = 


(1.2) 


'■'■■■'.5?  f-  ■';  rp-' 
Here  p  is  the  density,  m  and  E  are  the  momentum  and  total  energy  per  umt  volume  respec- 
tively and  p  is  the  thermodynamic  pressure.  These  equations  represent  the  laws  of  conserva- 
tion of  mass,  momentimi  and  energy  respectively.  The  thermodynamic  variables  p,  p  and  E 
are  related  by  a  caloric  equation  of  state  p  =  p(p,  £).  For  the  case  of  a  polytropic  gas  this 
relation  is  given  by: 


P  =  (-y-i) 


E  - 


m 


2p 


where  -y  is  positive  constant  usually  satisfying  K-y^— . 

Much  progress  has  recently  been  made  in  adapting  a  front  tracking  method  to  the  calcu- 
lation of  solutions  of  system  (1.1)  which  contain  discontinuities.   In  this  method  a  one- 
dimensional  grid  is  placed  onto  the  discontinuity.  Points  on  the  tracked  front  are  propagated 
by  solving  one-dimensional  Riemann  problems  in  the  direction  normal  to  the  interface.  This 
step  provides  the  position  of  the  tracked  interface  for  the  next  time  step.  Tangential  informa- 
tion is  ignored  during  the  normal  propagation  phase,  so  this  step  is  followed  by  a  update  of 
the  states  on  the  new  interface  based  on  the  component  of  (1.1)  tangent  to  the  interface.   The 
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positions  of  the  new  and  old  fronts  together  with  their  assigned  states  are  used  as  boundary 
value  data  for  the  solution  of  the  states  off  the  front.  A  detailed  description  of  this  method 
can  be  found  in  [1].  . 

The  discontinuities  supported  by  the  Euler  equations  (1.2)  are  of  two  types,  shocks  and 
contact  discontinuities.  If  combustion  is  considered,  a  third  type  of  discontinuity,  a  combus- 
tion wave,  may  also  occur. 

In  this  paper  we  will  report  on  recent  progress  which  has  been  made  by  the  authors  and 
co-workers  in  modeling  solutions  in  two  space  dimensions  of  the  Euler  equations  for  detona- 
tion waves,  surface  instabilities  for  non-combustion  interactions  and  the  numerical  solution  of 


certain  two-dimensional  Riemann  problems.'  In  addition  recent  theoretical  work  concerning  a 
general  theory  of  elementary  waves  and  Riemann  solutions  for  systems  of  hyperbolic  conser- 
vation  laws  will  also  be  discussed.  '        I 

2.   Multi-dimensional  Riemann  Problems  and  Elementary  Waves 


While  the  theory  ot  hyperbolic  conservation  laws  in  one  space  dimension  is  highly 
developed,  tiie  corresponding  theory  for  two  or  more  space  dimensions  is  not  so  well  under- 
stood.  Recent  work  has  been  devoted  to  the  development  of  some  basic  notions  which  can 
be  used  in  an  S  matrix  theory  for  systems  of  hyperbolic  conservation  laws  in  more  than  one 
space  dimension,  see  [2,3,4]. 

Let  A  =  be  the  Jacobian  matrix  of  F.   Equation  (1.1)  is  Sciid  to  be  hyperbolic  in  a 

da 

domain  D  Q  R"'*'^  ii  the  sxs  matrix  /fl  has  real  eigenvalues  X^ X,  for  all  (x,  r)  c  D 

and  for  all  vectors  |  .   If  the  eigenvalues  X^  are  all  distinct,  the  equation  is  said  to  be  strictly 
hyperbolic.   In  the  discussion  which  foUows,  hyperbolidty  is  assumed. 

An  S  matrix  theory  is  concerned  with  the  large  time  asymptotic  behavior  of  solutions  to 
systems  of  equations  for  which  system  (1.1)  is  a  first  approximation.   The  leading  order 
terms  of  these  large  time  solutions  are  governed  by  the  infinite  scaling  limit  of  the  original 
system  of  equations.   This  scaling  generally  eliminates  the  higher  order  effects,  yielding  the 


system  (11).  It  is  assumed  that  any  source  terms  in  the  original  equation  are  of  bounded 

'  ji  B  n:  r:..- 

extend.  Under  scaling  these  source  terms  will  in  general  survive  and  go  into  a  multiple  of  a 
delta  function  at  the  origin. 

An  S  matrix  is  the  product  of  two  wave  operators,  the  outgoing  operator  W*"  which 
gives  the  large  positive  time  asymptotics  and  the  incoming  wave  operator  W~  which  gives  the 
large  negative  time  asymptotics.   Attention  will  be  focused  on  the  outgoing  wave  operator 
W*".  The  domain  of  W^  usually  taken  to  be  the  range  ofW~.  However,  because  of  the 
occurrence  of  shocks  in  solutions  of  (1.1),  this  equation,  when  supplemented  by  the  necessary 
entropy  condition  to  separate  physical  from  nonphysical  waves,  is  not  reversible.  Thus  W~  is 
not  well  defined.   As  a  substitute  the  domain  of  W*"  will.be  restricted  to  scale  invariant  func- 
tions.  Thus  we  consider  solutions  to  the  initial  value  problem  for  system  (1.1)  whose  initial 
data  is  constant  on  rays  through  the  origin.  In  son;te  cases  it  will  also  be  desirable  to  impose 
regularity  conditions  on  the  initial  data  as  well.  In  one  Sjpace  dimension  this  is  the  well 
known  Riemann  problem  and  the  problem  of  solving  (1.1)  :with  scale  invariant  data  will  be 
referred  to  as  a  multi-dimensional  Riemann  problem. 

The  notion  of  dimensionality  in  a  Riemann  problem  is  actually  best  described  in  terms 
of  a  co-dimension.  A  Riemann  problem  of  co-dimension  j  is  defined  as  the  Cauchy  problem 
for  a  system  of  conservation  laws  in  d  space  dimensions  in  which  the  data  is  scale  invariant  in 
j  dimensions  and  independent  in  the  remaining  d  -  j  dimensions. 

The  restriction  to  scale  invariant  data  and  the  fact  that  equation  (1.1)  is  itself  scale 
invariant  implies  that  a  solution  to  a  Riemann  problem  should  be  self-similar,  that  is,  a  func- 
tion of  — .  This  implies  that  a  Riemann  solution  n  of  (1.1)  will  satisfy 

-  -Vo  -t-  VF  =  0.  (2.1) 

Such  a  solution  a  is  completely  determined  by  its  values  in  the  hyperplane  r  =  1  and  by  res- 
tricting our  attention  to  this  hyperplane  time  can  be  eliminated  from  the  equation.   Therefore 
a  Riemaim  solution  of  (1.1)  has  in  general  one  less  degree  of  freedom  than  a  general  solu- 
tion. 
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General  solutions  for  Riemann  problems  are  known  in  a  few  specisil  cases.   If  the 
number  of  space  dimensions  is  one,  the  system  is  strictly  hyperbolic  and  each  eigenvalue  is 
either  genuinely  non-linear  or  linearly  degenerate,  then  the  classical  paper  of  Lax   [5] 
describes  the  solution  of  a  Riemann  problem  with  a  small  discontinuity  as  consisting  of 
shocks,  centered  rarefactions  and  contact  discontinuities.   If  the  equation  is  scalar  and  in  two 
space  dimensions,  solutions  are  known  in  the  case  where  F  is  of  the  form  /V,  where  /  is  a 
scalar  valued  function  with  at  most  one  inflection  point  and  i.'  is  a  constant  vector  [6,7]. 

A  central  aspect  of  a  scattering  theory  is  that  a  source  decomposes  into  some  number  of 
localized  coherent  waves,  which  then  separate  and  propagate  away  from  each  other.   These 
local  disturbances  are  called  d- dimensional  elementary  waves  if  the  equation  (1.1)  is  in  d 
space  dimensions.   When  d  equals  one  these  elementary  waves  include  shocks,  contact 
diiicontinmties  and  rarefaction  waves.   For  d>\  erementary  waves  are  defined  by  the  interac- 
tion of  lower  dimensional  waves,  for  example  when  two  shocks  or  a  shock  and  a  contact 
discontinuity  collide.   In  many  cases,  such  as  the  interaction  of  two  shock  waves,  these  waves 
will  move  with  a  definite  velocity.   Assuming  that" the  original  equation  (1.1)  is  invariant 
under  Galilean  boosts  one  can  then  make  a'fcranslation  to  a  reference  frsime  in  which  the 
wave  is  at  rest,  thus  eliminating  one  more  degree  of  freedom  from  the  equation. 

Understanding  the  structure  of  these  elementary  waves  is  crucial  in  describing  the  solu- 
tion of  a  Riemann  problem.   In  some  special  cases  this  structure  is  known.   As  mentioned 
above,  in  the  case  of  one  space  dimension,  hyperbolicity  and  a  suitable  type  of  convexity  for 
the  flux  function,  a  theory  of  Lax  describes  these  elementary  wave  as  consisting  of  shocks, 
contact  discontinuities  and  rarefaction  waves  which  correspond  to  the  eigenvalues  of  the  dif- 
ferential of  the  flux  function.   If  the  equation  is  scalar,  then  a  general  theory  of  Oleinik  [8] 
describes  these  waves  in  terms  of  the  convex  envelope  of  the  associated  flux  function.   Other 
special  cases  in  which  the  structure  of  elementary  waves  is  known  include  one  space  dimen- 
sion polytropic  gas  dynamics  [9],  two  space  dimension  polytropic  gas  dynamics  [10],  adsorp- 
tion with  a  Langmuir  isotherm  [11],  and  water  and  polymer  displacement  of  oil  without 
adsorption  [12, 13]. 
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3.  Two  Dimensional  Detonation  Fi»9Dte>  r: 

The  method  of  front  tracking  has  been  applied  to  detonation  waves  in  two-dimensional 
gas  dynamics  [14].   For  problems  which  exhibit  cylindrical  symmetry  comparisons  can  be 
made  between  the  two-dimensional  model  and  a  one-dimensional  model  which  exploits  the 
symmetry  and  the  agreement  between  these  two  methods  is  good.   One  finds  that  as  the  mesh 
of  the  computational  grid  is  refined,  the  two-diraenr.io7ial  model  converges  linearly  to  the 
solution  given  by  the  model  based  on  cylindrical  symmetry. 

Only  a  polytropic  equation  of  state  is  conskdered,  so  the  energy  term  £  in  (1.2) 
becomes:  :-f.^  sbubd  saixrgi. 

•.::■--:.:;  3f{T    J^'     . 

Here  q  is  the  the  energy  released  by  tlie  chemical  reaction  that  occurs  across  the  detonation 
fro?n.  The  Chapman-Jouguet  model  of  detonation  tseused;  In  this  model  it  is  assumed  that 
the  reaction  takes  place  instantaneously  and  that  the  i^eaction  zone  is  infinitely  thin. 

If  state  0  is  the  unbumed  gas  and  state  1  is  theljurTSed  gas,  the  states  on  the  two  sides  of 
the  detonation  are  related  by: 

M2  =   ELZH^  (3.13) 

where  M  is  the  mass  flux  across  the  front,  and  the  Hugoniot  relation: 

In  the  case  of  a  Chapman-Jouguet  (CJ)  detonation,  the  detonation  wave  moves  at  the 
local  sound  speed  with  respect  to  tl^.e  gas  behind  it  [15],  and  the  behind  state  is  completely 
determined  by  the  ahead  state  and  equations  (3.1).  When  the  combustion  front  is  a  strong 
detonation  one  additional  parameter  is  necessary  in  order  to  specify  the  behind  state  from  the 
ahead  state. 


A  series  of  both  strong  and  CI  detonation  runs  using  grid  sizes  of  5  by  5,  10  by  10,  20 
by  20,  40  by  40  and  80  by  80  have  been  made.  The  contact  discontinuity  behind  the  detona- 
tion front  and  the  detonation  front  itself  were  tracked.   Several  of  these  runs  were  initially 
cylindrically  symmetric  and  in  these  cases  comparisons  were  made  with  a  one-dimensional 
computation  using  the  random  choice  method  with  1500  points  in  the  radial  direction. 

Figs.  3.1  -  3.5  present  the  results  of  a  cylindrically  symmetric  computation  in  which  the 
initial  pressure  ratio  across  the  front  is  100.   The  initial  dsnsity  is  uniform  and  the  gas 
releases  92.65%  of  its  internal  energy  upon  t-ombustion.   Figs.  3.1  and  3.2  show  the  positions 
of  the  contact  (inner  qucirter  circle)  and  the  detonation  (outer  quarter  circle)  at  the  beginning 
and  end  of  the  run  respectively.   Other  figures  include  comparisons  of  pressure  profiles  and 

detonation  wave  speed  (see  Figs  3.3  and  3.4).  The  detonation  wave  speed  error  when  calcu- 

■  1'' 

lated  with  respect  to  the  one-dimensional  code  is  Ifes's  than  0.5%.   Fig.  3.5  shows  the  conver- 
gence of  the  front  tracking  code  to  the  one-dimensional  code  under  mesh  refinement. 

In  addition  to  cylindrically  symmetricruns  the  front  tracking  code  has  2i]so  been  applied 
to  problems  in  which  the  initial  interface  is  elUpticJ.   If  the  initial  states  are  the  same  as  the 
ones  described  above,  hot  spots  are  produced  behind  the  front  in  regions  of  small  curvature 
and  cold  spots  in  corresponding  regions  of  large  curvature.   The  initial  lengths  of  the  major 
and  minor  axes  are  .3  and  .15  for  the  detonation  wave  and  .29  and  .145  for  the  contact.   Fig. 
3.6  shows  pressure  contours  and  the  waves  just  before  the  detonation  wave  breaks  through 
the  boundary  on  a  30  by  30  grid.   The  pressure  is  higher  behind  the  flatter  portion  of  the 
detonation  wave  than  behind  the  rounder  portion  of  the  wave. 


4.  Supersonic  Interface  Instabilities 

Interface  fingering  instabilities  arise  in  a  wide  variety  of  physical  contexts:  inertial  laser 
fusion,  plasma  fusion,  instabilities  of  layers  in  stars,  the  instability  of  laser  accelerated  foils, 
and  astrophysical  jets.   We  have  examined  the  compressible  Rayleigh-Taylor  instability  of  a 
superscnit  accelerated  contact  discontinuity  between  two  gases.  The  computed  solutions 
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exhibit  a  complicated  set  of  nonlinear  waves  comprised  of  spike  and  bubble  bow  shocks,  ter- 
minal shocks  within  the  spike  and  bubble,  ICelvin-Helmholtz  roll-up  of  the  spike  tip,  and  con- 
tact surface  waves.  Detailed  analysis  is  given  in  Ref.   [16].   We  have  also  studied  the  propa- 
gation of  a  supersonic  slab  jet  in  order  to  compare  and  contrast  the  jet  wave  structure  with 
that  of  the  supersonic  accelerated  sirrface. 

A  compressible  gas  interface  which  is  accelerated  by  a  shock  (the  Meshkov-Richtmyer 
instability  [17, 18])  is  Rayleigh-Taylor  unstable.  K  the  interface  is  accelerated  by  a  gravita- 
tional field,  then  the  interface  is  unstable  when  the  light  fluid  pushes  the  heavy.  The  impor- 
tant features  of  this  instability  can  be  modeled  by  impsrtijig  an  initial  kinetic  energy  to  the 
contact  discontinuity,  which  subsequently  u>  zUowe^to  ?dvsct  fieely.  We  assume  that  the 
problem  is  periodic  in  x  with  reflecting  boundaries  a('tBe  lop  and  bottom  of  the  computa- 
tional region.  .,:  .3,  rj.^jg  qj  oc; - 

The  problem  can  be  parametrized  in  dimensionless  units  by  the  initial  Mach  number  of 
the  tip  of  the  spike  with  respect  to  the  heavy  gas  and  by  theinitial  density  ratio  pj/p^  (b 
denotes  the  gas  below  the  contact,  a  the  gas  above).  Thfe^fflensional  scales  are  set  by  the 
initial  ambient  pressure,  perturbation  wavelength,  and  initial  ambient  density  of  the  heavy 
gas.  The  polytropic  gas  constant  -y  was  set  equal  to  1.4. 

An  interesting  set  of  wave  structures  emerges  from  this  study.  Figure  4.1  portrays  the 
evolution  of  a  Mach  2.8  density  ratio  2  accelerated  surface  at  r  =  0.4.  The  flow  is  initially 
supersonic  in  both  gases.  The  bow  shocks  in  the  lower  gas  have  interacted  to  form  a  single 
shock,  while  the  spike  bow  shock  has  interacted  and  joined  with  its  periodic  neighbors.  The 
spike  exhibits  the  characteristic  Rayleigh-Taylor  roll-up,  and  the  contact  shape  indicates  the 
presence  of  small-scale  surface  instabilities. 

Just  inside  of  the  advancing  spike  a  "terminal"  shock  wave  is  formed.  The  contact  is 
advancing  more  slowly  than  the  heavy  gas  inside  of  the  spike.   A  shock  wave,  preceded  by  a 
rarefaction  wave,  is  formed  as  the  advancing  heavy  gas  is  slowed  down  to  the  contact  velo- 
city.  A  similar  terminal  wave  is  formed  in  the  light  gas  inside  of  the  advancing  bubble. 

This  shock  preceded  by  a  rarefaction  wave  pattern  can  be  clearly  seen  in  the  density 
cross  section  plots  in  both  the  supersonic  accelerated  surface  run  (Figure  4.1)  and  the 
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supersonic  jet  run  (Figure  4.2).  Note  that  while  the  jet  terminal  shock  propigates  with  the 
head  of  tiie  jet  beam,  the  accelerated  surface  terminal  chc»cks  are  physical  transients  which 
decouple  from  the  late  evolution  of  the  contact  instability. 

The  compressible  Rayleigh-Taylor  results  differ  from  the  incompressible  case  chiefly  in 
the  formation  of  the  terminal  compression  waves  and  in  the  fact  that  the  spike  exhibits  less 
roll-up.   The  accelerated  surface  problem  differs  from  the  gravitational  instability  in  that  the 
spike  appears  to  attain  a  finite  growth  of  aspect  ratio  approjcimately  equal  to  2  for  cur  range 
of  parameters. 

A  resurgence  of  interest  in  supersonic  jtts  has  been  sparked  by  the  observation  of  astro- 
physical  jets  emanating  from  the  corca  of  adive.'gaiaxies.  aid  by  the  subsequent  success  of 
theoretical   [19]  and  computational  Mii'l}St&{20)i.  >...l^:,^. 

The  evolution  of  a  Mach  3,  density  ratio  10  slab  jet  at  r  =  0.4  is  presented  in  Figure  4.2 
[16].   The  jet  was  initialized  by  injecting  gas  at.a-specified  Mavh  number  into  an  a^abient  gas 
at  equal  pressure.   The  boundary  conditions,  are  ithrough-flow.   The  problem  is  parametrized 
by  the  Mach  number  of  the  jet  with  respect  to  the  jet  gas  and  the  density  ratio  of  jet  to 
ambient  gas.  -y  was  set  equal  to  5/3.   The  results  apply  to  jets  from  laboratory  to  astrophysi- 
cal  scales  since  the  problem  is  independent  of  length  scale. 

The  jet  bezim  in  our  80x120  grid  computation  is  5  grid  blocks  wide,  while  the  beam  is  20 
grid  blocks  wide  in  the  160x300  grid  computation  of  Norman,  Smarr,  and  Winkler  [20]. 

The  density  contour  and  cross-section  plots  in  Fig.  4.2  indicate  the  presence  of  a  bow 
wave  (the  flow  is  subsonic  in  the  ambient  gas)  and  of  a  terminal  shock  near  the  head  of  the 
jet  beam,  preceded  by  a  rarefaction  wave.  This  terminal  shock  system  may  explain  the 
observed  hot  spots  terminating  astrophysical  jets  [20].   The  contact  shape  displays  the  large 
scale  Kelvin-Helmholtz  roll-up  of  this  jet,  and  the  development  of  two-dimensional  pinch 
waves. 

The  fact  that  we  get  reasonable  results  with  a  beam  5  grid  blocks  across  illustrates  one 
of  the  advantages  of  the  front-tracking  method.   By  placing  additional  degrees  of  freedom 
around  the  tracked  contact,  the  method  is  able  to  resolve  the  solution  globally  with  fewer 
degrees  of  freedom  than  required  by  conventional  finite  difference  methods.   The  importance 
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of  this  feature  of  the  method  will  become  apparent  when  the  statistical  regime  of  multiple 
fingers  is  considered.  -?  : ■  .\ 


5.  Numerical  Implementation  of  Elementary  Waves 

Further  work  on  the  development  of  computer  code  for  modeling  elementary  waves  is 
in  progress.  Previous  papers  [10, 1]  have  reponed  the  numerical  implementation  in  the  front 
tracking  code  of  the  elementary  waves  known  as  regular  reflection  and  single  Mach  reflec- 
tion, this  section  will  discuss  the  case  of  shock  andicontact  discontinuity  interactions. 

The  simplest  model  of  a  shock  and  contact  diSi*)ntiiluity  interaction  consists  of  an 
incident  shock  wave  colliding  with  a  contact  d?s<ft)btihui^  separating  two  different  gases.  The 
local  result  of  this  interaction  is  a  com'iguradon  we  call  adiffraction  node.   A  diffraction 
node  consists  of  the  incident  shock  wave,  the  contact  dis'coritinuity  into  which  the  incident 
shock  collides,  a  reflected  wave  which  is  either  a  shoek  dr  a  centered  rarefaction  wave  and  a 
deflected  contact  surface  behind  the  incident  shock.  The  model  supposes  that  locally  all  of 
the  shock  or  contact  waves  can  be  assumed  to  be  straight,  and  that  the  solution  in  a  neighbor- 
hood of  the  node  is  piecewise  constant  except  for  the  possible  reflected  centered  rarefaction 
wave.   It  is  assumed  that  the  point  of  intersection  of  these  waves  moves  with  a  definite  velo- 
city and  thus  the  interaction  can  be  studied  in  a  frame  ol"  reference  in  which  the  node  is  at 
rest.  The  description  of  a  diffraction  node  then  consists  of  the  states  in  a  neighborhood  of 
the  node  together  with  the  angles  at  which  the  waves  at  the  node  intersect. 

In  a  dynamic  model  it  is  necessary  to  calculate  the  transformation  to  the  steady  frame  of 
the  node.   This  is  equivalent  to  finding  the  velocity  of  the  node  in  the  given  reference  frame. 
This  velocity  can  be  approximated  by  propagating  the  incident  shock  and  the  contact  discon- 
tinuity into  which  it  collides  for  one  time  step  ignoring  their  interaction.  The  intersection  of 
the  two  propagated  curves  is  then  used  as  the  updated  node  position  from  which  the  node 
velocity  can  be  calculated.   Once  the  node  velocity  is  known,  the  transformation  to  the  steady 
frame  of  the  node  is  performed.   K  one  assumes  that  the  data  in  front  of  the  incident  shock 
on  both  sides  the  contact  discontinuity  is  known,  and  the  strength  of  the  incident  shock  is 
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given,  then  the  configuration  around  tne  diffractiofl  ncde  in  the  s-.eady  frame  can  be  found  by 
the  intersection  of  shock  polars  in  the  pressure  turning  angle  space,  see  Henderson  [21]. 
Finally  the  configuration  is  translated  back  to  the  original  frame  of  reference. 

Figure  5.1  shows  the  result  of  the  interaction  of  a  planar  shock  wave  colliding  with  a 
sinusoidally  perturbed  contact  discontinuity.  The  incident  shock  is  in  air  (-y  =  1.402)  and  the 
contact  surface  separates  air  from  sulphur  hexafluoride  (-y  =  1.092).  This  interaction  is 
known  as  a  fast-slow  interaction  since  tiie  sound  speed  in  air  is  greater  than  that  in  sulphur 
hexafluoride.   The  initial  shock  has  a  pressui^  behind  to  pressure  ahead  ratio  of  100.   It  is 
interesting  to  note  the  extreme  proximity  of  tlie  transmitted  shock  and  the  deflected  contact 
discontinuity  near  the  node.   We  weri  ^^jjl-i  pltes'ed  with  ^Sl^  fiont  tracking  code's  ability  to 
resolve  a  configuration  with  such  close  cjrve.«.    The  rectang'ilar  mesh  used  for  this  nm  was 
20x20,  and  the  separation  between  the  transmitted  shpck  and  deflected  contact  is  on  the  order 
of  one  tenth  of  a  mesh  block  for  a  large  portion  of  the  computational  region.   Furthermore 
the  transmitted  wave  lies  on  the  sulphur  hexaflouride  side  of  the  contact  discontinuity  and  the 
value  of  gamma  for  this  gas  is  so  close  to  one  as  tp. make  the  resolution  of  waves  on 
moderate  sized  grids  difficult  for  most  finite  difference  methods.   We  suspect  that  without 
front  tracking  one  would  need  a  rectangular  raesh  more  than  ten  times  as  fine  in  each  linear 
dimension  in  order  to  resolve  both  the  deflected  contact  and  the  transmitted  wave. 


6.   Bifurcations  and  Wave  Interactions 

One  of  the  principal  difficulties  in  any  front  tracking  code  is  the  handling  of  bifurcations 
and  interactions  in  the  tracked  interface.   Examples  of  these  include  the  transition  from  regu- 
lar to  Mach  reflection  and  the  passing  of  waves  through  computational  boundaries.   Interac- 
tions of  both  of  these  types  have  been  either  fully  or  partially  implemented  in  our  front  track- 
ing code. 

Figure  6.1  shows  a  planar  shock  wave  incident  upon  a  ramp.   The  problem  is  initialized 
with  the  shock  normal  to  the  wall.   When  the  ramp  is  reached,  a  bifurcation  occurs,  in  this 
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case  to  a  regular  reflection.  At  a  later^time  the  point  of  regular  reflection  reaches  the  top  of 
the  ramp.  At  this  point  the  regular  reflssdoa  node  lifts  off  the  wall  producing  a  Macb  type 
reflection. 

In  figures  6.2  the  front  tracking  code  is  wsed  to  follow  the  development  of  a  compressi- 
ble Kelvin-Helmholtz  roll-up  [1,22].   Initially  two  gases  of  equal  pressure  and  temperature 
but  moving  in  opposite  directions  are  separated  by  a  slip  discontinuity.  This  slip  surface  is 
given  £in  initial  perturbation  which  causes  it  to  roU-up.  The  boundary  conditions  at  the  sides 
of  the  computational  rectangle  are  periodic.  As  the  surface  rolls  up  portions  of  the  surface 
cross  the  periodic  boundaries.   Any  section  of  the  surface  which  propagates  past  a  periodic 
boundary  is  disconnected  from  the  original  curve  and  reinstalled  periodically  shifted  to  the 
opposite  side  of  computational  rectangle.   A  linking  between  the  periodically  connected 
curves  is  maintained  *o  that  periodic  boundary  conditions  are  enforced.  The  visual  effect  is 
that  as  one  section  of  the  slip  surface  propagates  out  of  the  computational  window,  we  see 
the  corresponding  portion  of  the  periodic  neighbor  moving  into  our  pictxire. 

Not  only  is  the  interaction  of  curves  with  computational  boundaries  of  interest,  but  the 
interaction  of  nodes  as  well.  Fig.  6.3  shows  a  diffraction  node  propagating  past  a  computa- 
tionally passive  boundary.   This  problem  is  supersonic  and  the  signals  from  the  exterior  of 
the  right  hand  boundary  are  sufficiently  weak  that  their  influence  on  the  solution  can  be 
ignored.  Thus  the  problem  of  node  passing  through  such  a  boimdary  is  simply  a  matter  of 
identifying  those  curves  at  the  exiting  node  which  leave  and  those  which  remain.   Exiting 
curves  are  deleted  and  the  remaining  curves  are  separately  installed  on  the  boundary.   The 
main  difficulty  here  is  dealing  with  the  numerical  degeneracies  which  occur  because  very 
short  curves  are  produced  when  the  node  first  crosses  the  boundary. 
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Fig.  3.1 


The  initial  strong  detonation  wave  (outer  drcle)  and  the  contact  discontinuity  for  a 
cvUndrically  symmetric  computation.  The  initial  conditions  are  uniform  density,  zero 
velocity   and  a  circular  pressure  discontinuity  at  radius  .2,  with  ratio  inside  to  outside  of 
100    The  heat  released  upon  combustion  is  92.65%  of  the  internal  energy  of  the 
unbumed  gas.   The  initial  position  of  the  contact  is  radius  .195.   The  initial  state 
between  the  waves  is  that  behind  a  planar  detonation  wave  with  the  above  initial  data. 


Fig.  3.2 


The  detonation  wave  and  the  contact  at  the  time  step  analyzed  in  Figs.  3.3  and  3.4.  The 
detonation  wave  now  has  radius  approximately  .36. 


RADIUS 


Fig.  3.3 


A  plot  of  presure  vs.  radius  corresponding  to  Fig.  3.2  is  shown.   The  solid  curve  shows 
the  results  obtained  by  the  one-dimensional  random  choice  computation.   The  vertical 
lined  represent  the  range  of  pressure  values  in  the  two-dimensional  front  tracking  solu- 
tion at  a  fixed  radius  as  the  angle  varies  on  a  40  by  40  grid.    Thus,  the  vertical  lines 
show  the  anguJar  dependence  in  the  solution. 
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Fig.  3.4 


A  plot  of  detonation  speed  vs.  time  for  the  computation  on  a  40  by  40  grid  The  solid 
curve  shows  the  speed  of  the  detonation  wave  in  the  one-dimensional  cSculat^on  Te 
vertical  lines  represent  the  range  of  values  of  the  speed  of_the  detonation  in  the  two 

The  maximum  error  (max       'i'^^^)  is  less  than  0.5%. 


dimensional  calculation. 
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Fig.  3.5 

Convergence  of  the  front  and  interior  schemes.   The  pressure  errors  in  the  interior  and 
at  the  front  are  shown  for  NxN  grids  at  the  time  indicated  by  Fig.  3.2.   The  #  signs 
represent  the  interior  error,  where 


Interior  Error  =  100%  x 


Jo   /o   ^.nizHd^-J^y 


The  front  error  (error  bars)  gives  the  range  of  the  errors  at  the  front,  defined  as 

P    -P 

Front  Error  =  100%  x      -^       -^  , 

where  [P]  is  the  pressure  jump  at  the  front  in  the  one-diraensiona]  computation  at  the 
same  time.   The  asterisks  represent  the  error  of  the  average  pressure  behind  the  front, 
namely 


Front  Error(average  pressure)  -  100%  x 


P-  -P 

^d  averase  , 

[P] 
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Fig.  3.6 

Pressure  contours  are  shown  for  a  computation  of  an  elliptical  expanding  detonation  on 
a  30  by  30  grid.   Also  shown  are  the  detonation  wave  (D)  and  the  contact  (C). 
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Fig.  5.1 

The  interaction  of  a  shock  wave  with  a  contact  discontinuity  producing  reflected  and 
transmitted  shock  on  a  20  by  20  rectanguJar  grid. 
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Fig.  6.1 

A  shock  incident  upon  a  ramp.   Bifurcation  to  a  legulax  reflection  occurs  when  the 
shock  reaches  the  ramp.   When  the  regular  reflection  node  reaches  the  top  of  the  ramp 
a  bifurcation  to  a  Mach  type  reflection  occurs.  The  grid  here  is  30  by  30. 
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Fig.  6.2 

Compressible  Kelvin- Helmholtz  roll-up.  The  computation  is  on  a  40  by  40  grid.  When 
interior  curves  cross  the  periodic  boundaries  at  the  side  of  the  square,  they  are  periodi- 
cally reinserted  on  the  opposite  boundary. 
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Fig.  6.3 

Propagation  of  an  interior  node  past  a  computationally  passive  boundary  on  a  20  by  20 
grid.  Signals  from  the  outside  of  the  computational  domain  are  assumed  to  be  negligi- 
ble. 


Bukiet           •— "-'^"D      c.l 

Applications   of  front 
tracking   to   combus?Jon, . . 

NYU   DOE/ER/03077-266      c.l     _ 
Bukiet 
^  Applications   of   front              _ 

tracking  to   combustion... 

DATE   DUE 

BORROWERS    NAME 

Tliii  b<»(>k  mny  be  kt-pt 

FOURTEEN    DAYS 

A  fine  will  be  charfied  for  fQch  <iny  iht*  hrxik  is  tept  overtime. 

1        CAYUOnO  142 

